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ABSTRACT 


Subspace  methods  of  solving  spectral  estimation  and  direction  of  arrival  (DOA> 
problems  involve  finding  the  eigenvalues  and  eigenvectors  of  correlation  matrices.  Using 
the  Lanczos  algorithm  some  of  the  extreme  eigenvalues  and  eigenvectors  can  be  ap¬ 
proximated  without  requiring  the  entire  matrix  decomposition  theoretically  saving  many 
computations. 

This  thesis  develops  a  model  and  a  form  of  the  Lanczos  algorithm  to  solve  the  DOA 
problem.  The  relationship  of  the  number  of  eigenvectors  used  to  the  accuracy  of  the 
results  in  a  low  signal  to  noise  ratio  example  are  examined. 
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I.  INTRODUCTION 


Recently,  a  number  of  signal  processing  techniques  have  been  developed  that  in¬ 
volve  resolving  a  received  signal  into  so-called  signal  and  noise  subspaces.  The  ability 
to  perform  spectral  estimation  or  direction  of  arrival  from  the  determination  of  the  noise 
subspace  has  been  shown  in  the  works  of  Pisarenko,  Schmidt,  Kay,  Kailath,  and  others, 
The  methods  vary  in  the  manner  in  which  the  subspace  is  reached  and  how  the  resulting 
eigenvectors  are  calculated.  [References  1,2,  and  3J. 

To  achieve  better  results  (detection  at  lower  signal-to-noise  ratios,  better  resolution, 
finer  accuracy,  less  bias),  more  samples  are  required.  This  leads  to  larger,  mere  complex 
matrices  that  better  describe  the  signal  and  noise  subspaces.  Determination  of  the  sub- 
spaces  requires  an  eigendecomposition  of  an  estimated  correlation  matrix,  flic  general 
eigendecomposition  of  a  matrix  requires  computations  0[>r),  thus  the  larger  matrices 
require  many  more  operations.  Once  the  matrix  is  decomposed  into  its  eigenvalues  and 
eigenvectors  the  proper  set  of  eigenvectors  must  be  used  to  find  the  resulting  spectrum. 
Hence,  estimation  is  required  to  split  the  eigenvalues  into  signal  and  noise  subsets. 

The  difficulties  in  the  procedures  can  be  stated  with  two  questions. 

•  Where  is  the  threshold  between  signal  and  noise  eigenvalues  (and  thus  which,  or 

how  many  eigemectors  are  used); 

•  What  method  should  be  used  to  find  those  eigenvectors  in  a  timely  fashion? 

Proposed  here  is  a  method  which  will  accurately  estimate  some  of  the  eigem  allies 
and  eigenvectors  of  a  matrix  without  performing  the  entire  decomposition.  Computa¬ 
tional  savings  are  realized  when  only  a  partial  decomposition  is  required.  The 
eigenvectors  used  correspond  to  the  minimum  eigenvalues  of  the  matrix.  With  an 
over-specified  matrix  (dimension  much  larger  than  the  number  of  signals  present),  these 
minimum  eigenvalues  will  all  correspond  to  the  noise  subspace.  Mach  noise  eigenvector 
contains  all  the  information  to  find  the  actual  spectrum,  although  spurious  results  will 
also  be  apparent  (because  the  matrix  is  over  specified).  Several  methods  of  handling 
these  spurious  peaks  are  illustrated,  including  eigenvector  averaging,  spectral  averaging 
and  using  the  spectral  product. 

1  his  thesis  is  organized  in  five  chapters.  Chapter  II,  Direction  of  Arrival,  discusses 
classical  spectral  estimation  theory  and  how  it  applies  to  the  linear  array  problem 
(beamforming).  Subspace  methods  starting  with  Pisarenko  Harmonic  Decomposition 


and  proceeding  to  ML'SIC  and  I; SPRIT  are  discussed  in  detail.  Chapter  HI.  The 
Lanczos  Algorithm,  includes  all  the  basic  theory  required  to  describe  this 
eigendecomposition  method.  There  we  also  compare  several  methods  to  negate  the 
spurious  peaks  with  the  proposed  spectral  product  giving  the  best  results.  Results  of  the 
algorithm  for  numei  ous  cases  are  given  in  Chapter  IV.  The  last  chapter  summarizes  the 
results  and  advantages  of  this  Lanczos  subspace  method.  This  chapter  also  includes 
some  recommendations  for  future  work  and  lists  some  possible  applications. 


II.  DIRECTION  OF  ARRIVAL 

Direction  of  arrival  is  a  form  of  spectral  analysis,  performing  frequency  detection 
and  resolution  in  the  spatial  domain  vice  the  conventionally  considered  time  domain. 
The  signals  incident  on  an  array  are  analyzed,  and,  if  the  presumptions  of  the  analysis 
are  valid,  the  correct  bearings  to  the  sources  are  determined.  I  ormerlv  it  was  only  pos¬ 
sible  to  analyze  the  output  of  an  array  by  conventional  Fourier  techniques.  More  re¬ 
cently  numerous  methods  have  been  developed  which  enhance  the  ability  to  accurately 
determine  the  spectral  and  angular  resolution. 

This  chapter  summarizes  some  of  the  salient  points  of  spectral  estimation  be! ore 
discussing  classical  direction  of  arrival  array  processing  (Bartlett  beamforming). 
Projection-type  superresolution  subspace  methods  are  then  discussed,  starting  with 
Pisarenkos  Harmonic  Decomposition.  ML  SIC  and  FSPRi  I  are  discussed  in  detail  and 
several  other  subspace  methods  arc  mentioned. 

A.  SPLCTRAL  ESTIMATION 

Spectra!  estimation  is  the  term  used  to  describe  efforts  to  find  the  frequency  com¬ 
ponents  of  a  signal  sampled  in  time.  Two  conditions  that  arc  required  for  the  remainder 
of  this  thesis  are  that  the  processes  considered  are  wide  sense  stationary  (WSSi  and 
ergodic.  The  assumption  ol'WSS  means  that  the  process  has  finite  mean  and  that  tiie 
autocorrelation  is  a  function  of  the  distance,  or  lag.  between  two  samples  and  not  of  the 
position  of  the  samples  themselves.  Frgodicity  allows  time  averages  to  be  used  to  de¬ 
termine  ensemble  averages. 

The  autocorrelation  function  of  the  signal  \(t)  is 
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where  a in)  arc  the  individual  samples  of  the  signal.  When  only  a  finite  number  of  sam¬ 
ples.  ,\/,  are  taken,  the  above  definition  must  be  modified.  The  sample  autocorrelation 
function  is  defined  as 
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It  is  easily  shown  that  (Ref.  -4:  pp.  56-58] 

A„(  Ah  =  A„(  -A),  for  —  [M  -  1 )  <  A  <  0 

and  (3) 

Rxx(0)  >  Rxx(k),  for  all  A 


The  sample  autocorrelation  matrix  R,,  is 
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1  he  power  spectral  Jcnsin  <  PS  I ) »  is  a  measure  of  power  per  unit  frequency.  A  plot 
of  PSD  versus  frequency  will  show  the  relative  power  at  ail  the  frequencies  present.  It 
also  describes  the  properties  of  the  noise  in  the  signal,  i.e..  white  noise  will  have  a  fat 
PSD  showing  that  ail  frequencies  are  equally  represented  | Ref.  1:  pp.  53j.  !  he  PSD  is 

given  by 


•S' , , </:  =  ~ -  /"  R,r[m)c  ~/'r‘ 1 
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where  T  is  the  sampling  period. 

'1  he  periodograrn  method  of  estimating  the  true  PSD  is  one  of  the  earliest  and  most 
wide's  used  {Rei.  5:  pp.  5-Si.  I  he  periodograrn  is  defined  as  the  /-transform  of  the 
sample  auteemreh.lion  matiix  emulated  on  the  unit  siivle  i Re*.  pp.s'-'.y 
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It  may  also  be  found  by  directly  z-trar  forming  the  original  data  sequence  \(n) 


-prA'uU'tr  !)  where  A'(r) 


r.  =  I 


The  perioJogram  spectrum  is  found  by  substituting  r  =  c :,'T  . 
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Data  is  often  run  through  a  computationally  efficient  fast  Fourier  Transform  (FFT I  to 
find  the  perioJogram  spectrum. 

The  measures  of  effectiveness  of  a  spectral  estimator  are  based  on  comparisons  of 
resolution,  detectability,  bias  and  variance.  Resolution  relates  to  the  ability  of  the  esti¬ 
mator  to  separate  two  separate  spectral  lines  that  are  closely  spaced  in  frequency.  The 
capacity  to  locate  a  low  energy  signal  is  measured  in  detectibility.  Bias  is  the  tendency 
of  the  ooserved  peak  to  be  shifted  off  the  true  spectral  line.  Wiriance  is  a  measure  of  the 
propensity  to  keep  the  true  spectral  shape,  or  how  quickiy  the  PSD  lalb  oifaway  from 
the  true  peak.  Diflerent  spectral  estimators  are  sometimes  compared  in  terms  of 
robustness,  or  ability  to  function  weli  with  various  types  of  signals  and  noise. 

If  individual  signals  are  narrowband  the  resolution  criterion  is  said  to  be  achieved 
when  direct  onservution  of'  the  spectrum  leads  to  the  correct  determination  of  the  num¬ 
ber  of  signals.  Separate  peaks  are  not  required  to  deter:;  . me  that  two  signals  are  present 
[Ref.  6).  Resolution  is  inversely  proportional  to  the  length  of  the  data  samples.  With  f 
as  the  sampling  rate  and  T  the  sample  period,  or  time  between  samples.  7'=-yr.  the  fre¬ 
quent,  y  resolution  is  given  bv 


Ml 
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Thu*;,  with  the  periodogram.  better  resolution  can  only  be  achieved  with  longer  record 
lengths. 

The  above  description  is  based  on  no  windowing  (or  rectangular  windowing)  of  the 
data  samples.  The  windowing  in  the  time  domain  is  seen  as  convolution  in  the  frequence 
domain.  The  rectangular  windows,  for  example,  transform  into  sine  functions  in  the 
frequency  domain  resulting  in  high  sidelobe  distortion  known  as  leakage  in  the  frequency 
domain.  High  sidelobes  result  in  many  false  peaks,  making  it  difficult  to  decern  the  true 
spectral  peaks,  so  detectibility  suffers. 

Nenrectangular  windows  are  used  to  taper  the  data  to  minimize  the  discontinuities 
that  cause  the  high  sidelobes.  Common  windows  include  the  Bartlett  and  Hamming. 


The  lower  sidelobes  come  at  a  cost  of  resolution,  so  two  or  more  signals  dose  to  each 
other  in  frequency  may  merge  into  one.  Resolution  may  be  boosted  by  lengthening  the 


sample  time,  but  the  increased  record  length  may  s white  the  requirement  for 
stationarin .  More  can  be  found  on  the  subject  of  windows  in  References  "  and  S. 

Because  of  the  difficulty  m  meeting  the  requirements  for  both  detectibility  and  re¬ 
solution.  parametric  methods  have  been  developed  that  produce  increased  resolution 
with  short  data  lengths.  T  l,e  parametric  methods  are  based  on  determining  an  appro¬ 
priate  model  for  the  process  that  produced  the  data.  If  the  process  can  be  effectively 
modeled,  then  more  reasonable  assumptions  may  be  made  about  the  data's  behavior 
outside  of  the  window  ti  e..  these  data,  points  do  not  have  to  be  set  to  zero).  Once  the 
appropriate  model  is  chosen  to  represent  the  process,  the  parameters  arc  estimated  from 
the  available  data,  inserted  into  the  model,  and  the  power  spectral  density  expression  for 
the  respective  model  is  determined.  A  few  common  parametric  methods  are  summarized 
below.  [Ref  1 :  pp.  liK'-'us.  Ref 5:  pp.  172-174] 

Main  random  processes  that  are  encountered  in  practice  may  be  represented  by  the 
linear  difference  equation 
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where  u(n)  is  the  input  driving  sequence  and  x(n)  is  the  output  sequence.  The  a(k )  pa¬ 
rameters  are  the  autoregressive  coefficients  and  the  blk)  are  the  moving  average  cocfh- 
cients.  Equation  10  is  thus  known  as  the  autoregressive-moving  average  (ARMAi  model 
or  ARM  A ';>,«/  process  and  is  the  most  general  model  with  a  rational  transfer  function. 
The  power  spectral  density  that  results  from  the  ARM  A  model  exhibits  both  sharp 
peaks  and  deep  nulls.  If  the  autozcgtessive  parameters  of  equation  10  are  all  set  to  zero 
except  a'Oj  =  1,  the  resultant  model  is  the  moving  average  (MA)  process  of  order  q.  The 
MA  spectrum  will  have  deep  nn'ls,  but  relatively  broad  peaks.  With  the  b'kj  coeffi¬ 
cients  of  equation  10  set  to  zero  except  b‘0j  -  1.  the  autoregressive  (AR)  process  of  order 
p  results.  The  AR  spectrum  will  exhibit  sharp  peaks  but  will  have  broad  nulls.  (Ref  5: 
pp.  174-1  SI]  Each  of  the  models  will  exhibit  greater  accuracy  and  flexibility  as  the  order 
is  increased.  With  a  high  enough  order  the  AR  method  can  approximate  an  ARMA  or 
MA  process,  and.  likewise,  a  very  high  order  MA  model  can  be  used  to  approximate  an 
ARMA  cr  AR  process.  But  if  the  model  order  is  too  high  for  the  actual  process,  esti¬ 
mation  errors  will  lead  to  nonzero  coefficients  and  spurious  peaks  will  result  Thus 
proper  estimation  of  model  order  is  important  [Ref  1:  pp.  112-1 13.  pp.  l'kS-2<>"j. 

The  parameters  of  these  three  models  may  be  estimated  by  u-h.g  the  Yule-Walker 
equations  utilizing  the  autocorrelation  matrix  of  the  available  data  stream  [Re!  1:  pp. 
1 15-1  IS] 
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While  the  true  autocorrelation  function  is  impossible  to  determine,  the  maximum 
likclih»»ii  estimator  ( MI.  >  is  one  means  to  find  good  approximations  of  the  parameters 
for  tire  AR  model.  I  he  Ml.  method  uses  a  suitable  estimate  of  the  autocorrelation  or 
covariance  matrix  and  then  solves  [Ref  1:  pp.  ESb-ldO] 


(a  =  —  c 


for  the  parameters  a  where  C  is  the  symmetric  covariance  matrix  and  c  is  an 
autoeorrelated  vector. 

1  he  Burg  method  (maximum  entropv )  estimates  reflection  coelficients  from  the 
process  and  then  uses  the  Levinson  recursion  to  find  the  estimated  parameters  [Ref  1: 
pp.  22S-23 1  ]. 

Generally,  the  various  AR  spectral  estimators  except  the  autocorrelation  method  are 
unbiased  and  have  similar  variance.  The  covariance  and  Burg  methods  are  restricted  to 


ranges  that  keep  the  summation  within  the  available  data  and  do  not  assume  zero  pads 
outside  of  the  samples,  thus  taking  advantage  of  the  basis  which  led  to  the  creation  of 
the  parametric  methods  in  the  first  place.  [Ref  1:  pp.  240-252) 

B.  BEAMFORMING 

A  conventional  approach  to  the  direction  of  arrival  (DOA)  problem  is  through  the 
classical  beamforming  (Bartlett  )  technique.  Basically,  this  is  a  measure  of  coherency  of 
the  signals  arriving  at  an  array  of  sensors.  The  characteristics  of  an  array  are  described 
in  terms  of  array  gain,  directivity,  resolution,  beamwidth,  and  sidelobes.  These  are  based 
on  array  size  (number  of  sensors),  sensor  sensitivity,  intersensor  spacing,  and  post  re¬ 
ception  processing. 

Assuming  a  narrowband  point  source  in  the  far  field,  a  plane  wave  will  pass  through 
a  linear  array  in  a  specified  order,  where  the  magnitude  of  the  excitation  on  any  indi¬ 
vidual  sensor  will  be  related  to  the  phase  of  the  signal  at  the  instant  of  sampling.  1  he 
individual  sensors  will  see  different  instantaneous  magnitudes  based  on  this  phase  dif¬ 
ference  which  is  a  (unction  of  the  frequency  (or  wavelength),  the  DOA.  and  the  spacing 
between  sensors.  The  difference  in  plume  for  two  successive  sensors  for  a  linear  array 
with  equally  spaced  sensors  is 


A <i>  =  — -  sin  6  (15) 

/ 

where  d  is  the  distance  between  sensors.  /.  is  the  signal  wavelength,  and  0  is  the  angle 
between  the  wavefront  and  the  array  axis.  This  phase  difference  Ay<  is  also  known  as  the 
normalized  wavenumber.  A  [Ref.  4:  pp.  341-345).  Figure  1  illustrates  the  array • 
wavefront  interaction. 

The  output  of  a  simple  narrowband  delay-and-sum  beamformer  is 
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n.~  | 


where  .v„(f),  m=  1.  2. ...  ,  M  is  the  signal  at  the  with  sensor.  The  time  lag.  t„,  between 
two  adjacent  sensors  is  to  make  up  for  the  propagation  delay  caused  by  the  wavefront 
having  to  travel  the  extra  distance  d  sin  0.  One  can  adjust  the  magnitude  of  the  output 
to  plane  waves  arriving  from  a  particular  direction  0  by  introducing  appropriate  time 
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delays  at  each  sensor  prior  to  summing.  This  is  known  as  "steering  the  array"  or 
"steering  the  beam".  An  illustration  of  a  typical  beamformcr  arrangement  is  shown  in 
Figure  2. 


WAVEFRONT 


Figure  1.  Wavefront 

In  the  multiple  source  case,  especially  if  the  undcsired  signals  are  interfering  with  the 
detection  of  other  sources,  this  technique  may  be  modified  to  steer  nulls  instead  of 
beams  thus  minimizing  output  from  'jammers  '.  Nulls  may  be  directed  toward  a  single 
known  source  to  minimi/e  the  strength  of  the  signal  with  respect  to  that  source,  with  the 
output  being  analyzed  to  determine  what  other  sources  arc  present.  Another  imple¬ 
mentation  is  to  steer  the  nulls  to  minimize  the  total  output.  The  analysis  of  the  delays 
will  indicate  the  directions  of  multiple  taigcts.  jRcf.  4:  pp.  351-355.  Ref.  9  ] 

Beamforming  may  also  be  analyzed  in  the  frequency  domain  by  transforming 
equation  14  into  the  frequency  domain 


M 


m- 1 


(15) 


9 


In  vector  notation  we  can  write  e  =  wrx  where  w  are  the  weights  and  x  are  the  outputs 
of  M  sensors 
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The  steering  vector  s*  is  the  phase  relationship  between  the  angle  0  and  the  normalized 
wavenumber  k  given  in  equation  13 
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and  k  =  — —  sin  6  (17) 

/ 

I  ew-w  | 

It  can  be  shown  that  the  steering  vector  from  an  array  with  weights  w,  directed  toward 
an  arbitrary  direction  0,  can  be  expressed  in  terms  of  the  steering  vector  as  a  =  s.  [Ref. 
4:  PP.  34-344], 

Frequency  domain  beamforming  is  directly  analogous  to  spectral  estimation 
decribed  above.  Spatial  sampling  has  the  requirement  that  sensor  distances  d  must  be 
less  than  // 2  apart  to  prevent  "grating  lobes"  (or  aliasing i  corresponding  to  the  Nyquist 
rate  in  the  frequency  domain  of  /,ax  <J'j2  .  Longer  arrays,  containing  more  sensors, 
will  give  better  resolution  and  smoothing.  This  is  similar  to  frequency  resolution  being 
proportional  to  the  data  record  length  (a\/~  f  >•  1^°*-  PP-  341-345.  Ref.  1  < »:  pp. 
4-S,  Ref.  1 1  :  pp.  293-299] 

The  DO  A  is  actually  a  relative  comparison  of  observed  spatial  frequency  and  known 
signal  frequency.  The  spatial  frequency  is  greatest  on  endjlrc.  when  the  wavefront  is 
perpendicular  to  the  array  (0  =  90°  or  -ip).  Here  the  phase  difference-  between  adjacent 
sensors  is  at  a  maximum.  A  simultaneous  sampling  of  all  sensors  at  one  instant  in  time, 
or  snapshot.  will  indicate  the  maximum  spatial  frequency.  Observation  of  the  spatial 
wave  over  time  (with  a  known  temporal  sample  rate)  will  indicate  the  end  of  the  array 
where  the  source  is  located. 

On  broadside  iU  =  V  or  z »,  the  wavefront  excites  each  sensor  identically.  Spatial 
sampling  indicates  no  phase  difference  along  the  entire  array,  except  for  the  effects  of 
additive  noise,  resulting  in  the  computation  of  zero  spatial  liequency.  Unfortunately, 
linear  arrays  have  an  inherent  ambiguity  with  broadside  signals.  The  side  of  the  array 
at  which  the  source  is  located  cannot  be  determined  without  further  information.  Spatial 
frequency  is  illustrated  in  Figure  3. 

An  extra  requirement  in  the  standard  DOA  problem  is  a  priori  knowledge  of  in- 
coming  signal  frequencies.  Typically,  this  is  handled  via  a  bank  of  bandpass  filters  on 
the  output  of  the  sensors.  Data  streams  from  the  sensors  at  each  desired  center  (Ve- 


1 1 


0  5  »0  15  20  25  30 


ARRAY  SENSOR 


Figure  3.  Spatial  Frequency:  Two  signals  with  SNR  =  2dB, 

0,  =  1°  and  03  =  45°  for  two  snapshots  at  time  =  (a)  r„ ,  (/>)  Note 
the  variation  in  'DC  level'  as  the  snapshots  sample  the  nearly  broadside 
signal  at  different  phases. 
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quency  (/|,  f, ...)  are  processed  in  parallel,  resulting  in  spatial  frequencies  for  each 
temporal  frequency.  The  DOAs  are  calculated  by  comparing  these  spatial  frequencies 
with  the  center  frequencies  of  their  respective  filter  banks. 

Improvement  in  beamforming  may  be  seen  through  the  use  of  windows,  weighing 
each  sensor  output  by  the  appropriate  amount  to  narrow  the  beamwidth  or  lower 
sidelobes,  but,  as  discussed  earlier,  at  a  cost  of  lowering  overall  resolution. 

C.  SUBSPACE  METHODS 
1.  Introduction 

The  projection-type  subspace  method  utilizes  the  properties  of  the 
autocorrelation,  covariance,  or  modified  covariance  data  matrices  and  their 
eigenvalue  eigenvector  decomposition  into  signal  components  and  noise  components  in 
estimating  the  DOA.  Generally,  subspace  methods  use  an  assumed  property  of  the  data 
to  provide  good  resolution  if  the  data  fits  the  particular  assumptions,  even  for  extreme’}' 
short  data  sets.  If  the  data  tand  signals)  do  not  meet  the  assumptions,  the  results  tan 
be  quite  misleading.  The  assumptions  here  call  for  white  noise  and  a  signal  whose  esti¬ 
mated  autocorrelation  matrix  converges  to  the  true  autocorrelation  matrix  as  the  order 
is  increased. 

for  p  complex  sinusoids  in  additive  complex  white  noise  the  combined 
autocorrelation  function  of  the  signal  plus  noise  is  given  by 


R.prn)  -  )  Pc1' 
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We  define  R  as  the  signal  autocorrelation  matrix  of  rank  p.  R;;  =  n] I  of  lull  rank  M 
and  stive  the  autocorrelation  matrix  as 


= 


V  ,  H 
L  ‘  f- 


—  Rrr  +  R-- 


(19) 


where  I  is  an  M  by  M  identity  matrix.  R„  is  of  full  rank  M  due  to  the  noise  contribution, 
P  is  the  power  of  the  /th  complex  sinusoid,  a]  is  the  noise  variance,  and 
e,  =  Ql.  c'-'-'' . ]r  .  Unfortunately,  this  decomposition  is  not  possible  without 


absolute  knowledge  of  the  noise.  The  p  signal  vectors  e,  contain  the  frequency  infor¬ 
mation  and  are  related  to  the  eigenvectors  of  Rx>  by  v,  =  — =-  e,  for  i  <  p.  The 
eigendecomposition  of  R,x  is 


Rxx  =* 


(/,•  +  al)Vjxl!  + 


The  principal  eigenvectors  of  R„  arc  a  function  of  both  the  signal  and  noise.  If  no  signal 
is  present  the  autocorrelation  matrix  is  simply  a  diagonal  matrix  with  the  eigenvalues 
equal  to  the  variance  of  the  noise  [Ref.  1:  pp.  422-423]: 

a]  0 

0  o]  0  0 

0  •  • 

R,,=  (2i) 
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The  data  generated  from  taking  3/  samples  of  p  sinusoids  in  white  noise  can  be 
used  to  generate  an  autocorrelation  matrix  with  the  following  properties: 

•  The  theoretical  autocorrelation  matrix  will  be  composed  of  a  signal  autocorrelation 
matrix  and  a  noise  autocorrelation  matrix. 

•  The  signal  autocorrelation  matrix  is  not  full  rank  if  M>  p. 

•  The  p  principal  eigenvectors  of  the  signal  autocorrelation  matrix  may  be  used  to 
find  the  sinusoidal  frequencies. 

•  The  p  principal  eigenvectors  of  the  signal  autocorrelation  matrix  are  identical  to  the 
p  principal  eigenvectors  of  the  total  autocorrelation  matrix. 

The  matrix  will  have  a  minimum  eigenvalue  /.  =  a;.  [Ref.  1:  pp.  422-434] 

Furthermore,  the  noise  subspace  eigenvectors  are  orthogonal  to  the  signal 
eigenvectors,  or  any  linear  combination  of  signal  eigenvectors.  For  the  theoretical 
autocorrelation  matrix,  Rw.  all  eigenvalues  not  associated  with  the  signals  are  associated 
with  the  noise  and  arc  identical  in  magnitude  at  /.  =  a;  .  the  minimum  eigenvalue  of 
R.,;.  Thus, 


R.v/Vu  —  /T'.mv.’./  (22) 

The  zeros  of  this  minimum  eigenvector  polynomial 


v,0+  D--  ;  (23) 

j=o 

will  have  p  zeros  on  the  unit  circle  corresponding  to  signal  frequencies  [Ref.  A:  pp. 
335-337,  Ref.  5:  pp.  371-372]. 

For  the  simple  case  of  M  -  1  complex  sinusoids,  the  autocorrelation  matrix  R(, 
will  have  a  single  noise  subspace  eigenvector  vAf  with  its  associated  eigenvalue  /.  =  o;  . 
the  minimum  eigenvalue  of  R...  Thus,  the  resulting  zeros  correspond  exactly  to  the 
sinusoidal  frequencies. 

2.  ML'S  1C 

Mlitiple  Signal  Classification  (Ml  SIC)  is  a  form  of  a  noise  subspace  fre¬ 
quency  estimator,  based  on  noise  subspace  eigenvectors  with  uniform  weighting.  1  he 
ML  SIC  algorithm  finds  a  pseudo  spectral  estimate  from  [Ref.  5:  p.  3"3] 


ef  f) 


where  e  =  [1,  c  :" . c  :"  ]  '  The  advantage  of  this  method  is  in  its  generalized 

nature.  There  is  no  longer  is  a  requirement  for  uniform  spati.il  sampling.  Any  array 
geometry  will  provide  a  solution.  A  well-designed  array  will  eliminate  bearing  ambigui¬ 
ties  and  provide  unique  solutions  [Ref.  2:  pp.  ll)-2S[. 

R.O.  Schmidt  [Refs.  2,  12.  13]  has  shown  that  a  group  of  sensors  excited  by  a 
stationary  point  source  emitter  of  known  frequency  will  have  a  magnitude  and  phase 
relationship  (or  correspondence)  based  on  the  DOA  of  the  plane  wave.  This  corre¬ 
spondence  depends  on  the  array  geometry,  as  well  as  individual  sensors  characteristics 
(i.c..  directivity,  gain,  and  frequency  response),  and  may  not  be  unique  for  that  one  di¬ 
rection  of  arrival. 
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By  examining  the  signal  in  terms  of  an  M  dimensional  complex  field,  where  each 
sensor  provides  an  orthononnal  axis,  one  can  see  that  a  single  signal  will  result  in  one 
steering  vector.  This  steering  vector  describes  the  relationship  between  the  individual 
sensors  in  terms  of  phase  and  magnitude  differences,  thus  lor  any  unique  signal  fre¬ 
quency  and  direction  of  arrival  there  is  one  unique  steering  vector.  The  magnitude  of 
the  vector  will  vary  with  time,  but  its  direction  in  M  space  is  a  constant  determined  by 
the  amplitude  and  phase  relationship  of  the  sensors  for  that  particular  signal  as  illus¬ 
trated  in  Figure  4. 

The  theory  of  superposition  applies,  thus  with  two  signals  present  the  instanta¬ 
neous  received  steering  vector  is  a  linear  combination  of  the  individual  steering  vectors. 
The  time  varying  steering  vector  will  move  in  a  plane  that  is  spanned  by  the  individual 
steering  vectors.  Figure  5  shows  the  subspace  plane  spanned  by  two  steering  vectors. 
The  same  idea  can  be  expanded  to  a  case  ol'p  independent  signals  causing  the  steering 
v  ector  to  move  through  a  p  dimensional  subspace  (as  long  as  M  >  p). 

Unfortunately,  the  steering  vector  may  not  determine  the  actual  DOA  uniquely. 
In  the  example  of  a  one-dimensional  linear  array,  a  signal  gives  an  infinite  number  of 


Figure  4.  Steering  vector  for  3  sensors  and  I  signal 


Sensor  2 


Sensor  1 

Figure  5.  Signal  Suhspace  for  2  signals 

DOAs  that  lie  in  a  cone  that  is  formed  by  revolving  the  actual  DOA  about  the  axis  of 
the  array.  Thus  the  array  design  and  its  manifold  (expected  response)  will  play  a  large 
part  in  achieving  optimum  results.  The  array  manifold  describes  the  picdicted  response 
of  the  array  to  any  possible  signal  or  combination  of  signals.  I  he  manilold  may  be  es¬ 
timated  analytically  or  through  physical  calibration. 

Analytically  describing  an  array  is  a  complex  mathematical  procedure  for  all  but 
the  simplest  arrays  (i.e.,  equally  spaced  sensors  in  a  linear  array  or  a  3  element 
orthogonal  array).  It  also  assumes  that  absolute  knowledge  of  the  array  geometry  is 
available  --  an  assumption  which  can  lead  to  errors  if  the  differences  are  too  large. 

Calibration  with  test  sources  requires  a  known,  low  noise  environment  while  the 
test  sources  of  each  desired  frequency  arc  placed  in  a  finite  number  of  possible  locations. 
The  estimated  response  to  actual  signals  is  an  interpolation  of  these  responses.  Each  set 
of  calibration  parameters  requires  storage  in  memory;  this  results  in  overall  massive 
storage  requirements  for  a  comprehensive  grid. 

Several  parameters  such  as  the  number  of  signals  present,  the  directions  of  ar¬ 
rival  of  those  signals,  and  the  signal  strengths  can  be  determined  from  the  signal  sub- 
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space  information.  More  complex  models,  however,  can  be  developed  that  can 
determine  signal  frequency  and  polarization  as  described  in  References  2  and  13. 

We  will  now  describe  the  basic  steps  in  the  MUSIC  algorithm  for  the  DOA 
problem.  First,  we  sample  the  signals  and  compute  the  autocorrelation  matrix  of  the 
data  R.  Then,  we  determine  the  ordered  set  of  eigenvalues  and  eigenvectors  of  R.  In 
particular,  the  eigenvectors  associated  with  the  minimum  eigenvalues  must  be  accurate. 
In  the  theoretical  case  the  signal  eigenvalues  are  composed  of  signal  strength  { P  )  and 
noise  variance  (<?f)  and  the  nonprincipal  eigenvalues  will  all  be  identically  a]. 

We  now  determine  the  number  of  signals  by  eigenvalue  comparison.  A  simple 
counting  of  the  eigenvalues  greater  than  a]  will  give  the  number  of  signals  present  in  the 
theoretical  case.  However,  the  sample  autocorrelation  estimates  does  not  lead  to  such 
simple  results.  Small  power  level  signals  may  have  small  eigenvalues  (perhaps  smaller 
than  n]  >.  and  the  noise  eigenvalues  will  not  actually  be  identical  but  will  group  or  cluster 
near  an  approximation  of  m  .  Methods  of  solution  include  likelihood  ratio  tests  where 
the  gaps  between  the  eigenvalues  determine  threshold  placement  j Ref  2:  pp. 

Once  we  find  the  number  ol  signals  present  we  can  determine  the  signal  sub- 
space  by  the  span  of  the  first  p  eigenvectors 

\V  =  a  v-  .  \;]  (25) 

The  signal  mill  space  is  its  orthogonal  complement 

V,  ■  =  [v . v . .  ...  v .<]  ( 2('l 

We  now  find  the  intersection  of  the  signal  subspace  with  the  array  manifold.  1  he 
intersection  >  given  m  e^uotic n  2-t  'or  the  case  c:  tuc  uiiucrmh  spaced  linear  arrow  In 
actuahtv  the  intersection  is  estimated  with,  a  'best-lit"  method  used  to  determine  the 


optimum  res.: 


In  the  nonlinear  case  the  arro.v  manifold  is  much  more  diilict: 


t  to  rep- 


T wo  major  areas  ol  'difficulty  with  the  Ml'SK'  algorithm  are  the  calculation  and 
storage  ol  the  array  manifold  and  pcrlormmg  the  ei  gen  dec  ompo  si  t  i  o  n  of  the 
autocorreiation  matrix  that  resuits  from  very  large  arrays. 

3.  ESPRIT 

In  Reference  14.  Paulraj.  Roy.  and  Kailath  describe  Estimation  of  Signal 
Parameters  via  Rotational  Invariance  Tec!.;. nines  H  SPR1  IT.  a  <uhsp.:cC  method  which 
util./es  tv. o  identical  subarr.  vs  X  and  Y  A  known  distance  called  a  displacement  vector 


separates  the  two  parallel  suharrays.  hut  no  rotation  car.  he  presen'  loach  sensor  in  a 
matching  pair  (doublet  )  must  be  identical,  but  knowledge  of  individual  sensor  and  array 
response  characteristics  is  not  required. 

The  A’  elements  of  both  arrays  an.  sampled  simultaneously  with  the  signal  at 


each  sensor  beine  givai  bv 

c  ■"  _ r 

xfi)  =  V  sk{t)a{0k)  +  nm 
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P 


sm  0, 
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where  the  sampled  signal  at  each  sensor  in  a  doublet  differs  only  by  a  phase  term  and 
additive  noise.  With  p  signals  present,  v,  is  the  wavefront  at  the  reference  sensor  in  the 
X  array.  t>,  ;s  1 1  tc  1)0. \  relative  to  the  displacement  vector.  <;(#„)  is  the  response  cl'  die 
;th  seroor  in  a  subarray  relative  to  the  reference  sensor  in  that  array  for  a  signJ  from 
bearing  0..  J  is  the  magnitude  of  the  displacement  vectoi  and  n  the  noise  term,  lit  vector 
notation  th.e  signals  at  the  suKrr.tvs  are 


Mi)  -  Hsi/1  -  n, t  .0 
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The  vector  of  wavefronts  at  the  reference  sensor  in  array  X  is  represented  by  s. /).  All 
displacements  and  pi.ase  differences  are  based  on  this  sensor.  1  he  p  by  p  diagonal  ma¬ 
trix,  <1>.  contains  tiie  phase  delays  that  occur  with  each  set  of  doublets 

(b  =  Uiag[ <•  '.  <'  "J  (.'**) 

a  _  ; 

w  here  <;>.  =  sm  ,  as  shown  in  equation  I). 
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If  R.„  is  the  signal  autocorrelation  matrix,  the  subarray  autocorrelation  matrix 
is  given  by 

R„  =  HR,;Hr  *  a2vl  (31) 

The  cross  correlation  between  subarrays  is 

Rx,  =  HRp/(1)7H7  f32j 

where  H  is  the  direction  matrix  whose  columns  are  the  direction  vectors  for  the  p 
wavefronts 

h I0k)  =  lhidk).  tu(6k). ...  ,  hm(d,)Y  (33) 

Alter  some  manipulation  jRef.  15:  pp.  251-253],  we  obtain  the  matrix 

r  =  (R_W  —  )  -  ;  R  v, 

=  HR,.Hr- ,HR,.07H7  i  34) 

=  HRrr(l  —  ycI>r)M7' 

In  general,  there  will  be  />  eigenvalues  of  this  matrix.  But  when  y  =  c  the  it  It  row 

of  1 1  -  ;  'I>")  becomes  zero,  leaving  p- 1  eigenvalues.  Each  value  of  y  where  this  occurs  is 
called  a  generalized  eigenvalue  (GET.  Now,  the  DOAs  can  be  determined  by 

/■  ip-, 

0>.  -  arcsinl  — — -  )  (3n 

2.-.; 

Due  to  estimation  errors  in  the  calculation  of  the  correlation  matrices,  some 
error  will  be  induced.  Generally,  the  GT  s  will  be  moved  oil  the  unit  circle  and  out  from 
the  origin,  but  in  the  case  of  strong  signals,  they  will  still  be  easily  discernible.  It  should 
be  noted  that  poor  array  design  may  result  in  possible  ambiguities  (similar  to  Ml  SIC). 

Ad'.antages  to  note  o'er  the  Ml  SIC  algorithm  come  with  respect  to  the  array 
manifold.  With  ESPRIT,  no  manifold  is  required.  The  considerable  calibration.  elTorts 
and  storage  requirements  are  nonexistent.  However,  the  two  subarrays  must  be  identical 
in  ail  respects  and  must  be  positioned  exactly  parallel  to  each  oilier  [Ref  i-lj. 


4.  Other  Subspace  Methods 

Large  variance  effects  may  be  seen  in  the  above  methods  due  to  the  poor  reli¬ 
ability  of  the  estimation  of  the  eigenvectors  associated  with  the  minimum  eigenvalues 
of  the  estimated  autocorrelation  matrix  R„.  A  di'l'erent  method  of' spectral  estimation 
is  the  use  of  the  principal  components  where  only  the  eigenvectors  associated  with  the 
largest  eigenvalues  are  used.  Some  methods  have  tried  to  minimize  the  effect  of  noise 
by  using  R„  —  a;I  but  the  estimation  errors  in  both  R„  and  a ;  have  limited  their  suc¬ 
cess  [Ref.  1:  pp.  425].  Other  spatial  spectral  methods  include  the  parametric  methods 
such  as  AR  and  ML  [Ref.  1:  pp.  426-427). 

The  structured  matrix  approximation  of  Kumaresan  and  Shaw  [Ref.  16]  uses  K 
snapshots  in  time  of  an  .V  element  uniformly  spaced  linear  array  with  each  snapshot  in 
time  forming  a  column  of  a  data  matrix  X.  which  is  then  approximated  by  Xvf  in  the  least 
squares  sense.  The  bearings  information  is  then  calculated  using  the  properties  of  the 
Vandermonde  matrix. 

A  combination  of  frequency  domain  beamforming  and  autoregressive  modeling 
techniques  have  been  employed  in  Reference  17  to  estimate  the  direction  of  arm  a!  in  a 
multiple  source  localization  situation  using  a  planar  array. 

llalpeny  and  Childers  [Ref.  IS)  use  frequency-wavenumber  filters  to  break  the 
multiple  wave  case  down  to  a  succession  of  single  wave  problems. 

Reference  19  explains  an  algorithm  that  uses  non-noise  eigenvectors  from  a 
covariance  matrix  to  obtain  high  resolution  direction  of  arrival  for  narrow  band  sources. 

An  adaptive  beamforming  method  similar  to  a  minimum  energy  approach  is 
docribed  in  Reference  20.  The  eigenstructure  of  the  correlation  matrix  is  analyzed  and 
the  computations  are  modified  to  optimize  resolution  but  at  a  cost  of  arras  gain. 


III.  THE  LANCZOS  ALGORITHM 


Lanc7os  algorithms  provide  a  method  to  find  some  of  the  extreme  eigenvalues 
(smallest  or  largest)  and  their  associated  eigenvectors  from  large  matrices  with  fewer 
operations  than  required  in  an  entire  matrix  decomposition.  The  procedure  is  a 
tridiagonalization  of  the  original  matrix  based  on  an  iteration  developed  by  Cornelius 
Lanczos  [Ref.  21:  pp.  49-206J.  Once  the  matrix  is  in  a  tridiagonalized  form  the 
eigenvalues  can  be  easily  computed  using  a  symmetric  QR  algorithm  [Ref.  15:  pp. 
27S-279]  or  bisection  (with  or  without  Sturm  sequencing)  [Ref.  15:  pp.  305-30SJ.  The 
algorithm  takes  advantage  of  "minimized  iterations"  providing  quick  convergence  to  the 
final  tridiagonal  matrix  and  avoiding  accumulation  of  roundoff  error.  [Ref.  22) 

The  method  fell  from  favor  as  a  tridiagonalization  technique  with  the  advent  of  the 
Givens  and  Householder  transformations  later  on  in  the  1950s.  Givens  transformations 
[Ref.  15:  pp.  3S-47J  use  plane  rotations  (orthogonal  matrices)  to  zero  out  undesired  en¬ 
tries  in  the  matrix  undergoing  tridiagonalization.  The  Householder  methods  [Ref.  15: 
pp.  43 — 1~]  use  elementary  reflectors  to  accomplish  the  same  end.  Both  methods  are  in¬ 
herently  stable,  with  the  Householder  method  being  slightly  superior  in  both  speed  and 
accuracy.  Both  methods  outperformed  Lanczos  as  a  complete  tridiagonalization  proce¬ 
dure  in  the  general  case  where  all  eigenvalues  arc  required  [Ref.  23:  pp.  42-43]. 

The  real  power  of  the  Lanczos  method  however  lies  in  obtaining  the  extreme  values 
quickly.  The  entire  matrix  need  not  be  tridiagonalized  before  solutions  start  to  converge. 
Also,  if  the  original  matrix  is  sparse,  the  Lanczos  method  maintains  that  property,  sav  ¬ 
ing  even  more  computations.  Thus  for  large  matrices  (order  >  100)  the  Lanczos  method 
will  converge  on  extreme  eigenvalues  in  many  fewer  operations.  Recently  f  ufts  and 
Melissinos  [Ref.  24:  pp.  43-47 1  have  derived  a  variation  of  the  Lanczos  method  for  high 
resolution  spectral  estimation  and  showed  that  their  method  outperforms  both  the  sin¬ 
gular  value  decomposition  and  the  power  method  of  principal  eigenvector  and 
eigenvalue  computation.  Later  in  this  chapter,  it  will  be  shown  that  storage  require¬ 
ments  are  also  minimized. 

This  chapter  starts  with  an  explanation  of  the  classical  Lanczos  algorithm  as  devel¬ 
oped  by  Lanczos  and  modified  by  Paige  [Ref.  25J.  Then  we  will  discuss  die  unorthodox 
Lanczos  algorithm  of  Cullum  and  Willoughby  where  no  reorthogonalizaticn  is  per¬ 
formed  [Ref.  25).  finally,  the  algorithm  used  in  the  direction  of  arrival  problem  are 


discussed  in  detail.  Also,  we  present  some  results  of' the  algorithm  in  the  form  of  spectral 
estimation  of  multiple  tones. 

A.  CLASSICAL  LANCZOS  AND  ITS  VARIANTS 

The  original  algorithm  was  designed  as  a  means  to  directly  tridiagonalize  the  sym¬ 
metric  matrix  A.  The  development  of  Lanczos  algorithm  requires  the  knowledge  of 
Krylov  sequences  and  subspaces.  For  an  n  by  n  matrix  A  and  any  arbitrary  non-null 
n  by  1  vector  v,  the  Krylov  sequence  of  vectors  is  defined  as: 

vi+1  =  Ay  =  Av,  for  k=  1,  2,...,  n  (36) 

In  the  above  sequence  there  will  be  a  vector,  say  which  can  be  expressed  as  a 
linear  combination  of  a!!  the  preceding  vectors.  The  Krylov  matrix  of  rank  m  is  then 
given  by 

V„.  =  [v;  v2  v3  ...  v„.]  =  [v,  Av  .  A‘v; .  Aw-'v,j  (37 1 

and  the  Krylov  subspace  is  the  space  that  spurs  the  vectors  !y.  y,  ...  ,  vj, 

A'"'(A.v;)  -  >pan  [v,.  Av.,  A'v,.  ...  ,  A''~’v;j  (3N) 

The  columns  of  the  n  by  m  Krylov  matrix  V„  are  orthogonal.  Tiie  u  idiagonali/ation 
of  the  riv  en  matrix  A  is  then  obtained  as 


T  =  Y'aY„,  (39) 

where  T  is  an  m  by  m  tridiagonal  matrix.  Thus,  the  given  matrix  A  can  be 
tridiagonalized  provided  v.e  have  an  efficient  way  to  compute  the  orthogonal  matrix  Y„ 
or  to  compute  the  elements  of  matrix  T  by  performing  the  decomposition  of  equation 
39  directly  as  proposed  by  Lanczos  (Ref.  22.  Ref.  1 5J. 

The  most  direct  way  of  performing  the  tridiagonalization  assumes  that  T  =  V’AV 
where  V  =  fv,  y  ...  v„J  .  Note  that  A  is  a  symmetric  nbyn  matrix  and  V  is  an 
n  by  m  orthogonal  matrix  constructed  from  the  Krylov  sequence  of  vectors.  The  basic- 
recursion  for  a  real  n  by  n  matrb:  A  starts  with  a  randomly  generated  unit  Lanczos  vec¬ 
tor  v,.  Define  scalar  /;,  =  0  and  v,  =  0.  Define  Lanczos  vectors  v  and  elements  y.  and 
for  /  =  1,2 . M. 


/>.-;'.+ 1  =  Ay  -  3-,-V/  - 


(41) 


T 

a.  i  =  v.  Av; 

Pi+i  =  v(-i-i^vi  1^-i 


This  is  referred  to  as  the  basic  Lanczos  iteration  or  recursion.  The  actual  Lanczos  ma¬ 
trix  T„  j  =  1,2, ...  ,  M,  is  a  symmetric  tridiagonal  matrix  with  a, ,  1  <  i  <  j,  on  the  di¬ 
agonal,  and  ,  1  <  i  <  {j  —  1),  above  and  below  the  diagonal. 
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(43) 


Tach  new  Lanczos  vector  v.^  is  obtained  from  orthogonali/ing  the  vector  At\  with 
v,  and  v..,  .  The  elements  a,  and  /? are  the  scalar  coefficients  that  make  up  tlic  I.anc/os 
matrix.  The  basic  Lanczos  recursion  may  be  condensed  into  matrix  notation  by  defining 
V„  =  (v;.  v,. ...  ,  v.„).  Then  from  equations  40.  41.  and  42 


AV„.  -  V„.T„.  4- 


(44) 


where  e„  is  the  coordinate  vector  whose  mth  component  is  1  while  the  rest  of  the  com¬ 
ponents  arc  0.  T„  is  the  m  by  m  Lanczos  matrix,  the  diagonal  entries  arc 

T„l  A.A)  —  ak  for  1  <  k  <  m.  (45) 


and  the  entries  along  the  two  sub-diagonals  on  either  side  of  the  principal  diagonal  are 
Tm (/<./<  — r-  1 )  =  T„.( k  4  1  ,A)  ==  /q.+  .  lor  1  k  <  i?i  —  1  (4(0 

Note  that  A  is  never  modified  (unlike  in  Givens  and  Householder  transformations),  thus 
advantage  may  be  taken  of  any  existing  sparsity.  Also,  the  only  storage  requirements 
arc  for  the  three  Lanczos  vectors  (n  dimension),  the  Lanczos  matrix  T .  and  the  original 
matrix  A. 

We  summarize  the  actual  steps: 

•  Use  the  basic  recursion  of  equations  40.  41.  and  42  to  transform  the  symmetric 
matrix  A  into  a  family  of  symmetiic  tridiagonal  matrices  T  .  j  =  1.2 . M. 
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•  For  in  <  M  find  the  eigenvalues  of  T„,  n  (also  known  as  the  Ritz  values  of  T„>. 

•  Use  some  or  all  of  these  eigenvalues,  p,  as  approximations  of  the  eigenvalues  of 
A,  /  . 

•  For  each  eigenvalue  n  Find  a  unit  eigenvector  u  so  that  T„u  =  p u  . 

The  Ritz  vector  y  is  the  approximation  of  the  eigenvector  of  A.  It  is  found  from  map¬ 
ping  the  eigenvector  u  of  the  Lanczos  matrix. 

y  S  VMu  (47) 

So  the  eigendecomposition  of  T„  finally  results  in  the  Ritz  pair  (p,  y)  which  approxi¬ 
mates  the  eigenvalues  and  eigenvectors  of  A.  [Ref.  23:  pp.  32-33.,  Ref.  15:  pp.  322-325.] 

Unfortunately,  the  effects  of  finite  precision  arithmetic  prevent  the  theoretical 
Lanczos  algorithm  from  working.  The  eigenvalues  of  A  and  the  eigenvalues  of  T_  no 
longer  converge.  Roundoff  errors  are  partially  to  blame,  but  the  dominant  elTect  is  from 
the  loss  of  orthogonality  of  the  Lanczos  vectors  v  .  From  equation  40  it  can  be  seen 
that  a  small  will  have  great  effect  on  v_,.  Paige  showed  that  the  algorithm  runs 
within  allowable  error  constraints  until  it  starts  to  converge  on  the  first  eigenvalue.  At 
this  point  l'i  .  becomes  small  and  the  Lanczos  vectors  lose  orthogonality.  The  loss  of 
orthogonality  is  not  random,  however.  It  always  occurs  in  the  direction  of  the  Ritz 
vector  y. 

This  trouble  can  be  dealt  with  through  reorthogonali/ation.  But  questions  that  we 
must  answer  include: 

•  How  much  reorthogonali/ation  is  required? 

•  Where  should  it  be  performed? 

•  Reorthogonali/e  with  respect  to  what? 

Complete  reorihi'pniializiuion  is  the  first  choice,  inserting  a  Householder  matrix 
computation  into  the  Lanczos  algorithm.  This  enforces  orthogonality  among  all  the 
Lanczos  vectors  and  is  effective  at  keeping  the  system  stable.  But  the  extra  computa¬ 
tions  it  requires  negate  any  advantage  of  the  Lanczos  algorithm,  making  the  number  of 
computations  on  the  same  order  as  a  complete  decomposition.  Numerous  vectors  have 
to  be  stored  and  compared  requiring  many  swaps  into  and  out  of  storage.  [Ref.  15:  pp. 


Selective  reorthogonaliiation  is  clearly  more  efficient.  The  extra  computations  are 
performed  only  if  orthogonality  checks  indicate  loss  is  imminent.  Paige  shows  that  the 
loss  of  orthogonality  occurs  only  when  the  algorithm  converges  on  a  Rit/  pair.  Thus, 
instead  of  reorthogonalizing  against  every  other  Lanczos  vector,  using  the  less  numerous 
Ritz  vectors  will  be  sufficient.  Storage  is  therefore  minimized.  Only  when  all  eigenvalues 
of  the  matrix  are  required  does  this  method  become  computationally  more  expensive 
than  other  techniques.  [Ref.  26] 

The  last  option  is  no  reonhogonalization.  The  explanation  above  would  seem  to  in¬ 
dicate  that  maintaining  orthogonality  is  required.  However  by  analyzing  the  causes  and 
effects  of  the  original  loss  of  orthogonalization  one  can  insert  corrections  into  the  sol¬ 
ution  to  give  valid  eigenvalues  and  eigenvectors.  This  is  the  procedure  that  will  be  ex¬ 
amined  in  the  next  section. 

One  other  property  will  be  mentioned.  Here  the  single  vector  Lanczos  recursion 
described  above  will  not  find  duplicate  eigenvalues  of  the  matrix  A.  Parlett's  proof  uses 
the  power  method  to  expand  v  to  compute  the  columns  of  the  Krylov  matrix  K  (v.  A). 
If  J(\)  is  the  smallest  invariant  subspace  of  5v  which  contains  v.  then  the  projection  of 
A  onto  J'{\)  has  simple  eigenvalues.  Also,  the  Krylov  subspaces  fill  up  J(v)  and  for 
some  /  <  n  we  have 

span1,  v}  c=  A'2(A.  v)  c  A'?(A.  v)  c  ...  c  A'U,  v)  =  A'/+1(A.  v)  =  J{\)  (4K) 

The  result  is  that  some  eigenvectors  of  A  may  be  missed,  and  any  repeated  eigenvalues 
will  not  be  detected.  (Ref.  27:  pp.  235-239] 

The  multiple  eigenvalue  problem  can  be  treated  by  using  a  Block  Lanczos  algorithm. 
Instead  of  obtaining  a  triuiugonal  matrix,  the  result  is  a  banded  block  matrix,  where  the 
diagonal  is  M„.  an  /  by  /  matrix,  and  the  above  and  below  the  principal  diagonal  are 
upper  triangular  matrices  B,  .  Iiach  block  must  be  dimensioned  /  >  p  .  where  />  is  the 
estimated  multiplicity  oi  the  desired  eigenvalue. 

M,  B[  0 

B;  I\I2  Bj  0 

0  Bj  •  •  • 

0  •  •  •  0 

o  .  .  b/ 

0  B  m 
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This  is  analogous  to  the  general  case  (Ref.  15:  pp.  337-339].  The  block  Lanczos  algo¬ 
rithm  is  briefly  reviewed  at  the  end  of  this  chapter. 

The  above  discussion  assumes  that  the  given  matrix  A  is  real  and  symmetric.  Besides 
the  algorithms  summarized  in  this  section  for  a  real  and  symmetric  matrix  case,  there 
are  other  general  Lanczos  algorithms  proposed  in  the  literature  [Ref.  23)  that  are  suitable 
for  Hermitian  matrices,  non-symmetric  matrices,  and  for  rectangular  matrices. 

B.  IMPLEMENTATION 

The  Lanczos  phenomenon  states  that  a  few,  not  all,  of  the  eigenvalues  of  a  very 
large  matrix  A  can  be  computed  using  the  Lanczos  recursion  with  no 
reorthogonalization.  But  to  find  most  of  the  eigenvalues,  the  Lanczos  matrix,  T„,  will 
grow  in  size  approaching  that  of  A.  causing  the  loss  of  orthogonality  of  the  Lanczos 
vectors.  The  loss  ot  orthogonality  results  in  many  spurious  eigenvalues,  as  well  as  extra 
copies  of  good  eigenvalues.  In  any  case,  a  test  is  required  to  confirm  either: 

•  a  "found"  eigenvalue  is  good,  or 

•  the  eigenvalue  that  appears  is  spurious. 

Golub  and  Van  Loan.  Parlett.  and  Paige  [Refs.  15.  27.  and  28)  describe  procedures 
that  look  at  the  eigenvalues  for  each  T„  as  m  is  stepped  up  in  size.  All  the  eigenvalues 
of  T...  are  computed  at  each  step.  The  good  eigenvalues  will  repeat  at  each  larger  T^. 
while  the  spurious  eigenvalues  jump  around.  If  an  eigenvalue  does  not  match  at  con- 
secuti'.e  T.’s  it  may  be  considered  spurious  and  thrown  out.  If  a  good  eigenvalue  is 
mistakenx  deleted  (due  to  numerical  roundoif),  it  can  be  counted  on  to  reappear  in  the 
next  step. 

Cullum  and  Wiiloughby  [Ref.  23[  take  a  different  tack  by  developing  a  test  to  find 
and  eliminate  bad  eigenvalues  and  retain  the  rest.  The  advantage  here  is  in  utilizing  the 
machine's  tolerance  to  drop  bad  eigenvalues,  while  not  discarding  good  eigenvalues  that 
base  yet  to  comerge.  As  a  result  a  larger  spectrum  is  available  sooner,  even  though  it 
may  only  be  a  rough  estimate  of  where  the  eigenvalues  will  finally  converge. 

In  practice,  parts  of  the  Lanczos  recursion  (equations  41  and  42)  arc  replaced  by 

A  =  VfAv,  -  /!>,_,)  (50) 


->  *7 


and 


Computation  of  the  element  a,  is  a  modified  Gram-Schmidt  orthogonalization  proce¬ 
dure.  The  new  is  equivalent  to  the  of  equation  42  but  now  it  directly  controls 
the  size  of  the  Lanczos  vector. 

In  what  follows  we  describe  two  Lanczos  algorithms,  namely  the  single  vector 
Lanczos  algorithm  which  is  modified  and  analyzed  by  Paige  [Ref.  28|  and  the  block 
Lanczos  algorithm  described  by  Cullum  and  Willoughby  [Ref.  23].  Both  algorithms  have 
been  considered  for  the  estimation  of  the  directions-of-arrival  of  multiple  targets  in  noisy 
environments  in  this  thesis. 

1.  Single  Vector  Lanczos  Algorithm 

The  first  procedure  to  be  described  is  the  Paige's  single  vector  Lanczos  algo¬ 
rithm  [Ref.  2S]  for  real  symmetric  matrices.  The  single  vector  Lanczos  procedure  is  one 
of  the  most  straightforward  implementations  of  the  theory.  This  procedure  will  find 
some  or  many  of  the  eigenvalues  and  eigenvectors  of  a  real  symmetric  matrix  A  such  that 
Ax  = /x.  It  will  not  detect  repeated  eigenvalues.  However,  it  may  be  noted  that  for 
many  problems  of  interest  in  practice  we  do  not  have  strictly  multiple  eigenvalues.  Tor 
example,  in  the  direction-of-arriva!  problems  the  smallest  eigenvalues  of  the 
autocorrelation  matrix  corresponding  to  the  noise  associated  with  the  target  signals  are 
spread  over  a  small  range  rather  than  coinciding  on  the  same  value  [Ref.  2]. 

No  rcorthogonalization  is  performed  as  part  of  the  single  sector  Lanczos  al¬ 
gorithm.  As  mentioned  earlier,  the  Lanczos  vectors  begin  to  lose  their  orthogonality 
when  we  seek  to  estimate  all  or  most  of  the  eigenvalues  of  the  real  symmetric  matrix  A. 
Lor  the  application  under  consideration,  however,  wc  are  generally  interested  in  only  a 
few  of  the  minimum  eigenvalues  and  the  corresponding  eigenvectors.  It  is  mainly  for  this 
reason  that  wc  have  not  attempted  the  complete  or  selective  rcorthogonalization  of  the 
Lanczos  vectors  in  this  study. 

Now  we  shall  outline  the  basic  steps  involved  in  the  single  sector  Lanczos  al¬ 
gorithm.  T  his  is  based  on  the  recursion  described  by  equations  40,  50  and  51.  Based  on 
these  equations  Paige  [Ref.  28]  presented  four  different  single  vector  algorithms.  We 
have  adapted  one  of  in  this  study.  The  complete  eigenvalue  eigenvector  problem  actually 
consists  of  three  parts:  (a)  obtaining  the  tridiagonal  matrix  T_,  from  the  given  symmetric 
and  real  matrix  A  using  Paige's  recursion,  (b)  determining  the  smallest  eigenvalues  of 
T„  using  the  bisection  method  and  Sturm  sequencing,  and  (c)  estimating  the  corre¬ 
sponding  eigenvectors  by  computing  the  Ritz  vectors.  The  details  are  presented  in  the 


Step  1:  As  shown  in  equation  43  the  symmetric  tridiagonal  matrix  T  lias  entries 
a,  and  /?  along  its  principal,  and  the  adjacent  sub  and  super  diagonals,  respectively.  The 
following  recursive  expressions  are  then  used  to  compute  the  entries  of  the  tridiagonal 
matrix,  and  also  the  Lanczos  vectors  y,  [Ref.  2SJ: 

Initial  conditions:  v,  is  an  arbitrary  n  by  i  vector  such  that  |jv||2  =  1 

u,  =  Av, 

a0  =  Po  =  0 

for  j  -  1,2, ...  ,  m 

aJ  =  vIuj  <52) 

V  ,  =  U  :  —  ft:V* 

J  J  J  J 

Pj  =  i!«yll2 

')'+ 1  =  "jWj 

ll;+i  =  Av/+  i  ~  Pjvj 

where  w,  and  u  arc  some  intermediate  vectors.  The  vector  v,  is  obtained  by  filling  its 
entries  with  a  random  number  sequence  and  then  normalizing  it  with  respect  to  its 
Cuciidean  norm.  Now  T...  is  obtained  by  simply  filling  its  entries  as  in  equation  43.  Note 
that  in  <g  n  in  the  above.  One  quick  test  to  ensure  that  we  have  obtained  a  fairly  accurate 
estimate  ofT„  is  to  compute  the  product  v/Y  .  The  result  should  be  equal  to  <),  .  where 
C>„  is  the  Kroncckcr  delta  function. 

Step  2:  The  eigenvalues  of  the  tridiagonal  matrix  T„.  denoted  by  /<..  can  be 
computed  using  the  bisection  method  and  Sturm  sequencing.  Actually  one  could  obtain 
both  eigens allies  and  eigenvectors  ofT„  by  employing  such  methods  as  the  QR  algo¬ 
rithm.  However,  when  only  a  few  eigenvalues  are  required,  the  bisection  method  seems 
appropriate.  Tor  the  given  m  by  m  matrix  T,  we  define  the  characterstic  polynomial 

pt„(u)  =  det  (Tm  —  mI)  (33) 

which  can  be  recursively  computed  as  follows 

l\ ;(.u)  =  1 

p,(u)  =  a,  —  /( 

(34) 

p/u)  =  {y.j  -  u)p_.{u)  -  pj_ i/r_2(/i) 
for  j  =  2.  3.  ...  ,  m 
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The  roots  of  the  polynomial  pn{n)  are  the  required  eigenvalues.  For  our  appli¬ 
cation,  we  are  only  interested  in  a  small  range  of  eigenvalues  at  the  lowest  end  of  the 
eigenvalue  spectrum.  We  make  use  of  the  Sturm  sequencing  property  that  the 
eigenvalues  of  T._,  strictly  separate  those  ofTj  [Ref.  15:  pp.  305-307]  and  implement  the 
following  iteration: 

a  4-  b 
M  =  — j— 

b  =  n  if  pm(a)pjh)  <  0  (55) 

a  =  n  if  pja)pjb)  >  0 

and  we  repeat  the  above  as  long  as  ]  b  —  a  ]  >  e(  |  b  j  +  ]  a  '  ),  where  t  is  the  machine 
unit  round-off  error  and  [a./r]  is  the  range  of  our  required  eigenvalues.  Determining  the 
range  of  interest  in  our  application  may  require  some  a  priori  knowledge  about  the  signal 
to  noise  ratios  (SNR)  and  it  may  take  a  couple  of  iterations  to  do  this.  Some  alternatives 
to  the  iteration  given  in  equation  55  are  to  use  a  straightforward  polynomial  root  finder 
and  then  pick  the  roots  of  interest,  or  to  employ  the  well  known  L-D-U  factorization, 
both  of  which  may  not  be  computationally  efficient. 

Step  3:  There  arc  two  ways  to  obtain  the  eigenvectors  of  A  knowing  its 
eigenvalues.  /  .  Note  that  u.  are  the  estimates  of/...  In  the  first  method,  we  compute 
the  eigenvectors  of  .  denoted  by  t .  and  then  obtain  the  eigenvectors  of  A  given  by 

\  =  L„\.  (56) 

where  Lm  =  [v,  v,  ...  v_]  is  the  Lanc/os  matrix  which  ideally  is  the  same  as  the 
Krylov  matrix  of  equation  3~.  Note  that  i  n  .  t  )  e  T„;. 

The  second  method  involves  computing  the  Ritz  vectors  either  from  the 
Rayleigh  quotient  interation  or  by  the  orthogonal  iteration.  Here  we  assume  that  we 
have  good  estimates  of/  from  the  previous  step,  and  proceed  to  obtain  the  eigenvector 
x  by  minimizing  the  cost  function 

J  =  !’( A  -  /.jl)\j}2  (5") 

It  can  be  shown  that  a  simple  minimization  of  J  with  respect  to  x  yields  the  RaUcigh 
quotient  ofx 
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\TA X; 


T 

\  XV 


(5S) 


Therefore,  given  /  and  using  equation  5S  wc  can  formulate  the  Rayleigh  quo¬ 
tient  iteration  as  follows  [Ref.  15.  27] 


Initial  condition:  x0  is  an  arbitrary  vector  such  that  |jx0|!; 
for  k  =  0,  1,  2, ... 


1 


r(xd  = 


x[.Axk 


XX:, 


(591 


solve  (A  -  rfx^lDV;.^,  =  xk  for  \k± 


x,. 


Z+i 


i  I')  I 


where  y, . ,  is  some  intermediate  vector.  We  stop  the  iteration  when  converges  to 

a  constant  or  when  rix, )  one  of  the  known  eigenvalues.  At  each  iteration  step  we 
need  to  solve  an  tt  by  n  system  of  equations  in  this  method.  One  advantage  with  this 
method,  however,  is  that  it  converges  very  quickly.  Besides  the  above  iteration,  some 
other  methods  are  outlined  in  References  15  .  25  and  27.  We  remark  that  ifoniv  a  few 


eigenvalues  and  eigenvectors  <  say.  live i  are  required,  it  may  be  more  direct  to  use 
equation  5b. 

We  now  present  an  example  of' the  ability  of  the  single  vector  .  nc/os  algorithm 
to  estimate  the  direction  of  arrival,  or  to  find  spectral  lines  in  noise,  and  the  advantages 
in  extracting  more  than  one  eigenvalue  and  eigenvector  in  this  process.  We  consider  a 
signal  with  three  sinusoids  present  in  noise 


.v!  >/] 


V 

/Lt 


.  cos(  2 r/n)  4-  n(n I 


( ('< ) ! 


where  A  are  the  amplitudes  of  sinusoids,  f  are  the  normalized  spatial  frequencies 
(•i  iLf  <  <*.5  corresponding  to  .0  A  0  <  -p  )  .  and  n[n I  is  the  zero  mean  white  noise  with 


a  variance  ol  n:. 


We  have  computed  a  25  by  25  autocorrelation  matrix  of  .r(/7),  R„  ,  by  using  100 
data  samples.  We  have  used  the  covariance  method  for  this  purpose,  hence  R„  is  real 
and  symmetric.  The  eigenvectors  xt  of  R„  corresponding  to  the  lowest  eigenvalues  are 
computed  by  employing  the  single  vector  Lanczos  algorithm.  The  power  spectra!  density- 
estimates  are  computed  as  follows: 

2 

~ZT -  (61) 

V  ,V~'  :  ,  '«■' 

1=0 

where  x /(  are  the  elements  of  the  jth  eigenvector,  xy. 

Figure  6  and  Figure  7  show  the  power  spectral  density  (PSD)  estimates  of 
equation  60  with  an  SNR  of  10  dB  for  j  =  1  and  2,  respectively. 


Figure  6.  PSD  of  first  eigenvector 

Note  that  the  index  j  indicates  the  increasing  magnitude  of  the  eigenvalues.  Thus, 
(T. ,,  x,)  are  the  lowest  eigenvalue  and  the  corresponding  eigenvector.  In  both  figures  we 
have  the  peaks  at  the  correct  locations  (9°,  27°,  63°).  However,  they  both  have  spurious 
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smaller  peaks  at  different  locations.  We  can  observe  the  same  trend  for  the  first  five 
eigenvectors  as  shown  in  Figure  8,  where  the  spectral  estimates  are  over  laid  on  each 
other. 

Based  on  the  above  results  one  feels  that  we  could  employ  seme  kind  of  aver¬ 
aging  to  get  rid  of  the  spurious  peaks  and  improve  the  estimation  performance.  We  have 
implemented  two  such  methods:  the  eigenvector  averaging  and  the  spectral  averaging. 
Figure  9  shows  the  result  of  the  algebraic  averaging  of  the  first  5  eigenvectors,  and 
Figure  10  shows  the  result  of  the  algebraic  averaging  of  the  spectral  estimates  of  the 
same  eigenvectors.  As  seen  from  Figures  9  and  10,  eigenvector  averaging  yields  better 
results  than  spectral  averaging. 


Figure  7,  PSD  of  second  eigenvector 


Improved  results,  however,  were  obtained  by  using  what  is  called  spectral 
multiplication  which  is  obtained  by  taking  the  product  of  the  individual  spectra,  given 


s„  W  - 


where  J  is  a  predetermined  number  (J  <  m  <  n).  Figures  11  and  12  show  the  results  of 
spectral  multiplication  for  J  =  2  and  7=5,  respectively.  As  can  be  seen  in  these  fig¬ 
ures,  using  more  spectra  in  equation  62  greatly  improves  the  performance.  Also,  even 
for  J  =  2,  spectral  multiplication  outperforms  the  eigenvector  averaging  method. 


Figure  8.  Overlayed  PSDs  of  first  5  eigenvector 

In  the  remainder  of  the  thesis,  we  have  used  spectral  multiplication  in  prefer¬ 
ence  to  the  eigenvector  or  spectral  averaging.  Figure  13  shows  the  multiplication  of  5 
spectra  for  the  case  when  the  SNR  =  0  dB.  We  notice  a  sputious  peak  around 
6  =  74°.  More  spurious  peaks  are  observed  when  the  SNR  is  decreased  to  -3  dB  (see 
Figure  14)  and  Figure  15  shows  the  spectrum  for  the  SNR  but  we  have  used  the 
eigenvectors  6-10  in  this  case.  Improved  performance  is  obtained  as  shown  in 
Figure  16  (J  =  10)  and  Figure  17  (J  =  15)  by  using  more  and  more  eigenvectors  in  the 
spectral  multiplication. 

In  all  the  above  cases  we  always  observed  the  signal  spectral  peaks  at  the  right 
places.  The  spurious  peaks,  however,  did  not  appear  at  the  same  location  as  we  used  a 
different  eigenvector  to  compute  the  spectrum, 
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Figure  11.  Speclr.il  product  for  2  eigenvectors 


Figure  14.  Spectral  product  for  5  eigenvectors,  -5  dB:  Using  second  through  sixth 
eigenvectors 
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Figure  15.  Spectral  product  for  5  eigemectors,  -5  dB:  Using  sixth  through  tenth 
eigenvectors 


Figure  16.  Spectral  product  for  10  eigemectors 


Figure  17.  Spectral  product  for  15  eigenvectors,  -5  dB 


2.  Other  Methods 

The  single  vector  Lanczos  algorithm  will  not  determine  that  repeating 
eigenvalues  exist,  thus  it  cannot  find  the  corresponding  eigenvectors.  The  subspace  that 
results  has  an  incomplete  basis  as  it  is  described  only  by  the  eigenvectors  that  arc  com¬ 
puted.  The  solution  is  to  use  a  bloc!-,  method  that  is  analogous  to  the  single  vector 
Lanczos  algorithm.  As  we  mentioned  earlier,  the  block  form  of  the  Lanczos  algorithm 
does  find  eigenvectors  with  multiplicity  p  as  long  as  the  blocks  are  dimensioned  /  > /’. 
We  attempted  to  incorporate  the  Culluin  and  Willoughby  hybrid  block  Lanczos  algo¬ 
rithm  (Ref.  29  Chapter  8)  into  our  direction  of  arrival  model.  We  postulated  that  it 
would  be  desirable  to  compute  a  few  of  the  extreme  smallest  repeating  eigenvalues  and 
their  respective  eigenvectors.  However  we  were  never  able  to  get  the  program  to  reliably 
compute  good  eigenvalues  and  eigenvectors  for  the  autocorrelation  matrix.  This  has  not 
posed  a  problem  for  our  model  as  the  covariance  matrix  does  not  appear  to  have  re¬ 
peating  eigenvalues,  but  a  larger  order  matrix  may  indeed  include  duplicating  noise 
eigenvalues  and  require  an  algorithm  that  will  accurately  operate  with  that  perturbation. 
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The  algorithm  we  attempted  to  use  is  actually  a  hybrid  approach  to  finding  the 
eigenvalues  and  eigenvectors  of  a  real  symmetric  matrix  A.  For  insight  into  the  problem 
look  at  the  block  analogy  of  equations  4u,  41,  and  42.  Define  matrices 
B,  =  0  and  V0  =  0.  The  n  by  q  matrix  V,  has  columns  that  are  orthonormal  random 
vectors.  The  value  of  q  must  be  greater  than  or  equal  to  the  number  of  eigenvalues  to 
be  found. 


for  /  =  2,  3,  ...  ,  s 

r  (63) 

Vi+1B(+t  =  P;  =  AV,  —  V.M,  —  Vj_,  B,r 

Mi  =  V’/'(AVi-V|._IB  1)  (64) 

V(+1B/+1  -  P,  (65) 

The  matrix  B._ ,  is  a  modified  Gram-Schmidt  orthonorinali/ation  of  the  columns  of  P,. 
Also  note  that  the  matrix  M  correspond  to  the  a's  of  the  single  vector  Lanc/os. 

1  he  block  analogy  to  the  Krylov  subspace  approach  can  be  performed  with 

AlA.V.)  =  s/nnijv,.  AY;,  A*V>, ...  ,  A'-ivj  (66) 

The  blocks  V.  for  j  =  1.  2 . s  form  the  orthonormal  basis  of  the  Krylov  subspace. 

It  can  be  shown  that  for  a  symmetric  n  by  n  matrix  A  and  an  orthonormal 
n  by  q  starting  matrix  V..  that  the  block  recursion  equations  63,  64.  and  65  will  generate 

blocks  V,.  V. . V  where  qs  <  n.  It  is  these  blocks  that  form  an  orthonormal  basis 

of  the  subspace  A’(A.  V;)  .  In  much  the  same  way  as  the  single  sector  Lanc/os  algorithm 
generates  the  tridiagonal  Lanotos  matrix,  the  block  variant  generates  bl<wks.  but  these 
are  now  nontridiagonal.  At  the  end  of  each  iteration  the  Lanczos  matrix  is  of  the  form 
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A] 

1  b[ 

1 

m[am 

b2 

1  <*2  fa 

\  fa  a3 

fa 

m7am, 

1  fa 
\ 

•  •  • 

(67) 

1 

1 

•  •  • 

•  • 

fa 

1 

fa 

«x 

The  submatrix  MfAM  consists  of  the  reorthogonalized  terms  and  M.  is  the  portion  of 
the  first  block  that  is  not  generating  descendants.  Ritz  vectors  are  computed  on  every 
iteration  and  are  used  as  the  starting  blocks  for  the  next  iteration.  I'ach  block  is  required 
to  be  reorthogonalized  with  respect  to  the  all  the  vectors  in  first  block  which  is  not  being 
allowed  to  generate  descendants.  It  is  apparent  that  the  block  procedure  requires  a  great 
amount  of  storage  and  is  very  computationally  intensive. 
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IV.  RESULTS 


Using  the  Lanczos  algorithm  it  is  possible  to  line!  some  of  the  eigenvalues  and 
eigenvectors  of  a  matrix  without  going  through  an  entire  matrix  decomposition.  The 
smallest  eigenvalue  of  the  autocorrelation  matrix  and  its  corresponding  eigenvector  will 
have  the  required  spectral  information  to  determine  a  source’s  bearing  (direction  of  ar¬ 
rival)  from  an  array.  Multiplying  several  of  the  resultant  eigenvectors'  power  spectral 
densities  will  tend  to  reinforce  the  true  spectral  peak  and  zero  out  spurious  peaks  that 
do  not  occur  with  every  eigenvector. 

The  problem  with  finding  the  split  between  the  noise  and  signal  eigenvalues  disap¬ 
pears  as  only  a  few  of  the  smallest  eigenvalues  of  a  large  matrix  (in  relation  to  the 
number  of  sources)  are  used. 

A.  APPROACH 

The  received  signal  is  modeled  by  sinusoids  at  normalized  spatial  frequencies  pro¬ 
portional  to  their  bearings  from  cnJIire  (.0  =  0°,  .5  =  90°).  The  sum  of  these 
sinusoids  is  sampled  at  a  rate  based  on  the  interelement  spacing  of  2.  Thus  a  source 
at  endfire  is  sampled  at  the  Nyquist  rate  and  the  sample  rate  increases  as  the  bearing 
shifts  toward  the  array  broadside.  The  simulation  uses 


ss(n) 


A  cos(2rr f.n)  +  n{n) 


(6$) 


where  sstn)  is  the  instaneous  excitation  for  the  sensor  at  location  n.  A  is  the  amplitude 
of  each  of  the  T  signals.  ft  is  the  normalized  spatial  frequency  of  the  /th  source  (de¬ 
pendent  upon  bearing),  and  //(//)  is  white  gaussian  noise.  The  relationship  between  A 
and  the  noise  variance  a]  is  determined  by  the  desired  signal  to  noise  ratio  (SNR),  where 


SXR  = 


(69) 


The  experiment  consists  of  simulating  a  linear  array  with  equally  spaced  sensors  re¬ 
ceiving  signals  of  known  temporal  frequency  from  various  bearings.  One  possible 
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physical  implementation  would  place  a  bank  of  bandpass  filters  on  each  sensor  with  the 
outputs  from  each  similar  filter  tied  into  a  correlator.  Advantages  of  this  method  include 
*’.e  processing  gain  found  by  prefiltering  the  noise  and  simple  parallel  implementation 
with  separate  channels  for  each  frequency  band.  The  lowering  of  the  noise  bandwidth 
will  raise  the  SNR  at  the  correlator.  As  more  filters  are  used  (smaller  bandwidth)  the 
noise  power  decreases  and  the  SNR  is  increased.  The  algorithm  creates  an 
autocorrelation  matrix  with  the  output  of  the  correlator.  The  Lanczos  tridiagonalization 
and  eigendecomposition  provides  the  eigenvectors  that  are  estimates  of  the  spatial  PSD. 
The  PSD  corresponds  to  the  sources  directions  of  arrival.  A  possible  implementation  is 
shown  in  Figure  IS. 

B.  EXPERIMENT  SET  UP 

The  first  three  cases  show  the  effect  of  different  signal  strengths  on  the  ability  to 
accurately  determine  the  number  of  targets  and  the  bearing  resolution  for  various  di¬ 
rections  and  target  spacing.  In  each  of  these  cases,  the  number  of  sensors  is  100.  a  ma¬ 
trix  size  of  25  is  used  and  15  iterations  (the  number  of  eignevalucs  eigenvectors  found) 
are  performed.  Case  4  uses  three  5  dB  sources  at  18°,  36'  and  41.4°  to  illustrate  the  ef¬ 
fects  on  changing  the  number  of  sensors  (samples),  the  size  of  the  autocorrelation  ma¬ 
trix.  and  the  number  of  eigenvectors  used.  The  noise  is  randomly  generated  white 
gaussian  noise  with  a  standard  deviation  selected  to  provide  the  desired  SNR.  Each 
figure  shows,  uu  the  PSD  of  selected  eigenvectors  overlayed  and  plotted  versus  bearing 
and  (b)  the  product  of  selected  PSDs  of  those  eigenvectors. 

Case  1  is  with  all  sources  at  a  signal  strength  of  5  dB.  Figure  16  shows  results  from 
the  second  through  sixth  eigenvectors  of  a  single  source  at  IS0.  Note  that  some  of  the 
eigenvectors  have  individual  peaks  as  high  as  the  true  signal  peak,  but  only  at  the  true 
bearing  do  all  have  a  common  peak.  Figure  2V  illustrates  the  other  end  of  the  spectrum, 
at  ST.  Once  again  the  second  through  sixth  eigenvectors  are  overtaxed  to  show  that  the 
correct  bearing  is  consistently  displayed,  but  in  this  case  one  eigenvector  has  an  indi¬ 
vidual  peak  higher  than  the  true  signal  peak.  The  product  of  these  PSDs  provides  suf¬ 
ficient  resolution.  Figure  21  has  two  closely  spaced  sources  at  36°  and  3 S'.  Resolution 
is  achieved  but  the  PSD  product  of  the  second  through  sixth  eigenvectors  shows  a  spu¬ 
rious  peak  near  75°.  f  igure  22  is  from  three  sources  at  0°.  36°  and  SS.2°.  The  individual 
eigenvector  PSDs  clearly  show  the  excellent  performance  at  broadside. 

Case  2  lowers  the  signal  strength  of  all  sources  to  0  dB.  Figure  23  shows  results 
from  the  second  through  sixth  eigenvectors  of  a  single  source  at  IS0.  Many  more  peaks 


43 


Figure  18.  A  physical  implementation 

are  visible  in  the  PSD  product,  making  the  decision  of  how  many  targets  more  diflicult. 
Figure  24  shows  that  at  the  the  other  end  of  the  spectrum  (at  81°)  the  situation  is 


44 


Figure  20.  Case  1  5  dB,  1  target  at  81  (a)  overlay  of  PSDs  from  second 

through  sixth  eigenvectors,  (b)  product  of  the  above  PSDs 
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slightly  worse  (due  to  a  lower  sampling  rate).  The  overlays  of  the  second  through  sixth 
eigenvectors  show  that  the  correct  bearing  is  consistently  displayed,  but  in  this  case 
enough  spurious  peaks  reinforce  one  another,  resulting  in  the  PSD  product  that  has  not 
zeroed  out  the  bad  peaks.  Figure  25  illustrates  that  the  proper  choice  of  eigenvectors 
will  resolve  this  problem.  Here  the  PSDs  of  the  first  through  fifth  eigenvectors  are  used, 
giving  a  product  that  is  easier  to  determine  correctly.  Figure  26  shows  the  0  dB  case  for 
two  close  spaced  sources  at  36°  and  38°.  Resolution  is  achieved  but  the  PSD  product 
of  the  first  through  fifth  eigenvectors  shows  several  spurious  peaks,  including  the  same 
one  as  in  Case  1  near  75°.  Figure  27  is  the  three  source  example  at  0  dB.  The  individual 
eigenvector  PSDs  are  repeating  at  the  proper  bearings  but  the  performance  at  broadside 
is  resulting  in  the  product  at  the  other  bearings  actually  being  driven  down. 

A  signal  strength  of -5  dB  is  used  for  Case  3.  Figure  28  shows  results  from  the  first 
1<>  eigenvectors  of  a  single  source  at  1ST  F'sing  more  good  eigenvectors  increases  the 
likelihood  that  all  spurious  peaks  will  be  diminished.  Figure  29  illustrates  results  at  the 
ether  end  of  the  spectrum  (at  81°).  Resolution  is  looked  at  in  Figure  3<>.  At  -5  dB  the 
algorithm  cannot  separate  the  two  close  spaced  sources  at  36°  and  38s.  A  number  of 
spurious  pe..ks  arc  higher  than  the  bump  at  36°  making  it  impossible  to  accurately  de¬ 
termine  the  number  of  sources  as  well  as  both  locations.  Resolution  is  tried  again  in 
Figure  31  with  2  sources  at  and  40°  .  Using  five  eigenvector  PSD  products  produces 
good  results.  Figure  32  shows  3  sources  at  -5  dB.  Good  performance  is  seen  both,  in 
the  overlays  and  in  the  PSD  product. 

Case  4  starts  with  loo  sensors  and  a  25  by  25  covariance  matrix  shown  in 
Figure  33.  As  the  number  of  sensors  is  decreased  and  the  number  of  eigenvectors  used 
is  held  constant,  more  spurious  peaks  start  tu  occur  (Figure  34  and  Figure  35). 
Figure  36  shows  the  spectral  improvement  as  more  eigenvectors  are  used.  As  the 
number  of  sensors  is  decreased  to  40  (Figure  37)  and  seven  eigenvectors  arc  used,  the 
results  are  still  acceptable.  Using  only  3'  sensors  we  can  no  longer  see  resolve  between 
the  two  closely  spaced  sources.  With  a  sufficient  number  of  eigenvectors  no  spurious 
peaks  are  present  but  the  true  number  of  targets  is  nondeterminable  I  Figure  38  and 
I  igure  39). 
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Figure  24.  Case  2  0  (IB,  1  target  at  81  °:  (a)  overlay  of  PSDs  from  first  through 
fifth  eigenvectors,  (b)  product  of  the  second  through  sixth  eigenvector 
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Figure  25.  Case  2  0  dB.  1  target  at  81  °:  product  of  the  first  through  fifth 
eigenvector  I’SDs 
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Figure  26.  Case  2  U  dB.  2  targets  at  36  °  and  38  °:  (a)  overlay  of  PSDs  from  first 
through  third  eigenvectors,  (b)  product  of  the  first  through  fifth 
eigenvector  PSDs 
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Figure  31.  Case  3  -5  dB,  2  targets  at  36  0  ard  40  °:  (a)  overlay  of  PSDs  from  first 
through  sixth  eigenvectors,  (b)  product  of  the  first  5  eigenvector  PSDs 
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Figure  33.  Case  4  5  dB,  3  targets  at  18  36  0  and  41.4  100  sensors,  (a)  PSDs 

of  second  through  sixth  eigenvectors,  (b)  product  of  the  above  PSDs 
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of  second  through  sixth  eigenvectors,  (b)  product  of  the  above  PSDs 
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Power  Spectral  Density 


Figure  37.  Case  4  5  <IB,  3  targets  at  18  °,  36  0  and  41.4  °:  40  sensors,  20  by  20 
matrix,  (a)PSDs  of  first  through  seventh  eigenvectors,  (b)  product  of 
the  above  PSDs 
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Figure  38.  C  ase  4  5  dli,  3  targets  at  18  36  °  and  41.4  30  sensors,  20  by  20 

matrix,  (a)PSDs  of  first  through  fifth  eigenvectors,  (b)Product  of  the 
first  through  fourth  eigenvector  PSDs 
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Figure  39.  Case  4  5  dB,  3  targets  at  18  36  0  and  41.4  30  sensors,  20  by  20 

matrix,  (a)Product  oT  the  first  through  fifth  eigenvector  PSDs, 
(b)Product  of  the  first  through  twelfth  eigenvector  PSDs 


66 


V.  CONCLUSIONS  AND  RECOMMENDATIONS 


The  results  plotted  in  Chapter  IV  indicate  that  the  eigenvectors  found  using  the 
Lanczos  algorithm  are  sufficiently  accurate  to  determine  the  spectrum.  Although  no 
direct  comparisons  with  other  eigendccomposition  methods  are  performed,  the  theory 
indicates  that  many  fewer  operations  are  required.  We  handle  the  other  difficulty  of 
conventional  subspace  methods  by  using  only  a  few  of  the  eigenvectors  associated  with 
the  minimum  eigenvalues  of  the  autocorrelation  matrix.  No  estimation  of  the  noise 
subspace  dimension  is  required  or  performed. 

This  theory  may  be  applied  to  any  system  requiring  rapid  decomposition  of  the 
correlation  matrix.  Examples  include  phased  array  radar  and  passive  acoustic  arrays 
[Refs.  30.  31].  Reference  32  details  an  experimental  system  using  the  MUSIC  algorithm 
1'or  multiple  source  direction  finding. 

The  following  areas  are  recommended  for  future  study. 

•  Use  of  the  products  of  multiple  spectra  apparently  resulted  in  good  detection  at  low 
SNR.  More  research  in  this  area  to  determine  a  physical  interpretation  of  this 
method  is  required. 

•  Analysis  and  comparison  of  the  results  in  terms  of  computational  speed  and  accu¬ 
racy  with  other  eigendccomposition  methods  should  be  performed  to  find  the  true 
results. 

•  The  Lancvos  algorithm  developed  uses  no  reorthogonali/ation  nor  will  it  find  re¬ 
peating  eigenvalues.  Other  forms  of  the  Lancvos  algorithm  are  available.  Com¬ 
parisons  between  these  different  methods  to  determine  accuracy  and  speed  may 
lead  to  more  optimum  results. 

•  A  more  detailed  model  should  be  developed  that  will  simulate  an  array  with  a  hank 
of  bandpass  filters  to  better  forecast  results  of  a  physical  implementation  (as  seen 
in  Figure  Is  of  Chapter  IV). 

•  A  method  which  implements  the  algorithm  in  parallel  fashion  may  be  tried.  I  sing 
one  long  linear  array,  several  overlapping  subarrays  may  be  used  to  simultaneously 
create  several  autocorrelation  matrices.  The  algorithm  may  be  applied  to  these 
matrices  in  parallel.  It  is  predicted  that  the  greater  number  of  available 
eigenvectors  will  more  properly  describe  the  noise  subspace  and  therefore  more 
accurately  estimate  the  spectrum. 
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